#!/bin/csh
# timeseries_point.csh
####################################################################
# Plots timeseries of different fields at a given grid point
# It requires the timeseries_point.dat that can be generated with make_timeseries_point.dat
#
# - Set the simulation directory containing YACC's output 
#   in the file "simulation.dir"
# - Run the script using:
#   ./timeseries_point.csh time0 time1
#   where time0 and time1 are initial and final times of the series
#   If they are not given, the whole time series is plotted
#   By default, the timeseries of all fields marked with 1 are plotted
########################################################################

#-------------------------
# Possible time-series
#-------------------------
# Select timeseries to plot
set temperature = 1
set velocity = 1
set viscosity = 1
#-----------------------------------------------------------------------------------

#-----------------------------------------------------------------------------------
# Options:
#-----------------------------------------------------------------------------------
# dirout: directory of the YACC's output files (needs to be written in the 
#         'simulation.dir' file)
# dirplots: directory where to store the plots (default: $dirout/Plots/Timeseries)
# makejpg: make also jpg files (default: 1)
# dimensional_plot: if 1 the file scaling.dat must be in the simulation directory, otherwise 
#                   a non-dimensional plot is made (default: 0)
# jpgdensity: resolution of the converted jpg files (default: 150)
# deletepdf: delete pdf files (done only if makejpg=1, default: 1)
# dimx: plot width in cm (default: 20 cm)
# dimz: plot height in cm (default: 10 cm)
# X0: x-coordinate in cm of the bottom-left corner of the plot (default: 3 cm)
# Y0: Y-coordinate in cm of the bottom-left corner of the plot (default: 5 cm)
# filename: name of the file generated by make_timeseries_point.csh (default: timeseries_point.dat)
#
# Select options
set dimensional_plot = 1
set makejpg = 1
set jpgdensity = 200
set deletepdf = 1
set dimx = 20
set dimz = 10
set X0 = 3
set Y0 = 5
set filename = timeseries_point.dat
#-----------------------------------------------------------------------------------

source simulation.dir

if (! -d $dirout) then
  echo $dirout 'does not exist. Check the path stored in the variable "dirout" '
  exit
endif

set dirplots0 = {$dirout}Plots/
set dirplots = {$dirout}Plots/Timeseries/

#-------------------------------------------------------------------
# If the scaling file is not found, non-dimensional plots are made
#-------------------------------------------------------------------
if ($dimensional_plot == 1) then
  if (! -f {$dirout}scaling.dat) then
    echo {$dirout}scaling.dat 'for dimensional plot not found.'
    set dimensional_plot = 0
  endif
endif
  
if ($dimensional_plot == 1) then
  echo 'Dimensional plot'  
  source {$dirout}scaling.dat
else
  echo 'Non-dimensional plot'  
  set D_scale = 1.0
  set T_scale = 1.0
  set T_surf = 0.0
  set eta_ref = 1.0
  set k_ref = 1.0
  set alpha_ref = 1.0
  set time_scale = 1.0
  set V_scale = 1.0
  set eps_scale = 1.0
  set stress_scale = 1.0
  set heats_scale = 1.0
endif

#-------------------------------
# Select output files
#-------------------------------
set timeseries = ${dirout}${filename}

if (! -d $dirplots0) then
  mkdir $dirplots0
endif
if (! -d $dirplots) then
  mkdir $dirplots
endif
echo 'Plot grid-point time series for model:' $dirout

#----------------------------------
# Set gmtdefaults
# ---------------------------------
source setgmtdefaults

#-----------------------------------
# Set temporary ps output file
#------------------------------------
set output = timeseries_point.ps

#################
# Time-series 
#################
set time0 = ${1}
set time1 = ${2}
if (${1} == "" && ${2} == "") then
  set time0 = `awk 'NR > 1 {if(min==""){min="'"$time_scale"'"*$1}; if("'"$time_scale"'"*$1< min) {min="'"$time_scale"'"*$1}} END {printf "%12.5f", min}' $timeseries`
  set time1 = `awk 'NR > 1 {if(min==""){min=max="'"$time_scale"'"*$1}; if("'"$time_scale"'"*$1>max) {max="'"$time_scale"'"*$1}} END {printf "%12.5f", max}' $timeseries`
endif
set atime = `echo "($time1 - $time0)/5" | bc -l`

####################
# Temperature
####################
if ($temperature == 1) then
  set adaptive_scale = 0
  echo "Temperature point timeseries"

  set title = ""
  if($dimensional_plot == 1) then 
    set xlabel = "Time@~\040\133@~Ma@~\135@~"
    set zlabel = "Temperature@~\040\133@~K@~\135@~"
  else
    set xlabel = "Time"
    set zlabel = "Temperature"
  endif

  set pdfoutput = temperature_point_ts.pdf
  set jpgoutput = temperature_point_ts.jpg

  if ($adaptive_scale == 1) then
    set min = `awk '{if(min==""){min="'"$T_surf"'"+"'"$T_scale"'"*$7}; if("'"$T_surf"'"+"'"$T_scale"'"*$7< min) {min="'"$T_surf"'"+"'"$T_scale"'"*$7}} END {printf "%12.5f", min}' $timeseries`
    set max = `awk '{if(min==""){min=max="'"$T_surf"'"+"'"$T_scale"'"*$7}; if("'"$T_surf"'"+"'"$T_scale"'"*$7>max) {max="'"$T_surf"'"+"'"$T_scale"'"*$7}} END {printf "%12.5f", max}' $timeseries`
  else 
     set min = $T_surf
     set max = `echo "$T_surf + $T_scale" | bc -l`
  endif
  set delta = `echo "scale=1; ($max - $min)/5" | bc -l`

  psbasemap -X$X0 -Y$Y0 -R$time0/$time1/$min/$max/ -JX$dimx/$dimz -Ba{$atime}:${xlabel}:/a{$delta}:${zlabel}:WSen:.${title}: -K >! $output
  awk '{print "'"$time_scale"'"*$1, "'"$T_surf"'"+"'"$T_scale"'"*$7}' $timeseries | \
  psxy -R -JX -W5/0/0/0 -O >> $output

  ps2raster $output -Au -P -Tf -F$pdfoutput
  if ($makejpg == 1) then
    ps2raster $output -Au -P -Qt4 -Qg4 -E$jpgdensity -Tj -F$jpgoutput
  endif

endif

####################
# Velocity
####################
if ($velocity == 1) then
  set adaptive_scale = 1
  echo "Velocity point timeseries"

  set title = ""
  if($dimensional_plot == 1) then 
    set xlabel = "Time@~\040\133@~Ma@~\135@~"
    set zlabel = "Velocity@~\040\133@~cm/a@~\135@~"
  else
    set xlabel = "Time"
    set zlabel = "Velocity"
  endif

  set pdfoutput = velocity_point_ts.pdf
  set jpgoutput = velocity_point_ts.jpg

  if ($adaptive_scale == 1) then
    set min = `awk '{if(min==""){min="'"$V_scale"'"*$6}; if("'"$V_scale"'"*$6< min) {min="'"$V_scale"'"*$6}} END {printf "%12.5f", min}' $timeseries`
    set max = `awk '{if(min==""){min=max="'"$V_scale"'"*$6}; if("'"$V_scale"'"*$6>max) {max="'"$V_scale"'"*$6}} END {printf "%12.5f", max}' $timeseries`
  else 
     set min = 0
     set max = 10
  endif
  set delta = `echo "scale=1; ($max - $min)/5" | bc -l`

  psbasemap -X$X0 -Y$Y0 -R$time0/$time1/$min/$max/ -JX$dimx/$dimz -Ba{$atime}:${xlabel}:/a{$delta}:${zlabel}:WSen:.${title}: -K >! $output
  awk '{print "'"$time_scale"'"*$1, "'"$V_scale"'"*$6}' $timeseries | \
  psxy -R -JX -W5/0/0/0 -O >> $output

  ps2raster $output -Au -P -Tf -F$pdfoutput
  if ($makejpg == 1) then
    ps2raster $output -Au -P -Qt4 -Qg4 -E$jpgdensity -Tj -F$jpgoutput
  endif

endif

####################
# Viscosity
####################
if ($viscosity == 1) then
  set adaptive_scale = 1
  echo "Viscosity point timeseries"

  set title = ""
  if($dimensional_plot == 1) then 
    set xlabel = "Time@~\040\133@~Ma@~\135@~"
    set zlabel = "Viscosity@~\040\133@~Pa@~\040@~s@~\135@~"
  else
    set xlabel = "Time"
    set zlabel = "Viscosity"
  endif

  set pdfoutput = viscosity_point_ts.pdf
  set jpgoutput = viscosity_point_ts.jpg

  if ($adaptive_scale == 1) then
    set min = `awk '{if(min==""){min=log($8*"'"$eta_ref"'")/log(10)}; if(log($8*"'"$eta_ref"'")/log(10)< min) {min=log($8*"'"$eta_ref"'")/log(10)}} END {printf "%12.5f", min}' $timeseries`
    set max = `awk '{if(min==""){min=max=log($8*"'"$eta_ref"'")/log(10)}; if(log($8*"'"$eta_ref"'")/log(10)>max) {max=log($8*"'"$eta_ref"'")/log(10)}} END {printf "%12.5f", max}' $timeseries`
  else 
     set min = 1
     set max = 10
  endif
  set delta = `echo "scale=1; ($max - $min)/5" | bc -l`

  psbasemap -X$X0 -Y$Y0 -R$time0/$time1/$min/$max/ -JX$dimx/$dimz -Ba{$atime}:${xlabel}:/a{$delta}:${zlabel}:WSen:.${title}: -K >! $output
  awk '{print "'"$time_scale"'"*$1, log($8*"'"$eta_ref"'")/log(10)}' $timeseries | \
  psxy -R -JX -W5/0/0/0 -O >> $output

  ps2raster $output -Au -P -Tf -F$pdfoutput
  if ($makejpg == 1) then
    ps2raster $output -Au -P -Qt4 -Qg4 -E$jpgdensity -Tj -F$jpgoutput
  endif

endif

##############################
# Moving and cleaning
##############################

if ($makejpg == 1) then
  mv -f *.jpg $dirplots
endif

if ($makejpg == 1 && $deletepdf == 1) then
  rm -f *.pdf
else 
  mv -f *.pdf $dirplots
endif

rm -f *.ps 
